Superfluidity of fermion atoms loaded in a deep optical lattice: the existence of two 

rotonlike modes 
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We present theoretical calculations of collective modes of the one-band attractive Hubbard model 
which is widely used to study the s-wave superfluid phases of atomic Fermi gases of two-hyperfine 
states loaded in a deep optical lattice. To make our theory applicable for both superconductivity and 
superfluidity, we assume the more general t—U—J Hamiltonian. Using the functional differentiation 
we derive Schwinger-Dyson equations for the single-particle Green's functions. The method of 
Legendre transform is used to give a systematic derivation of the Bethe-Salpeter (BS) equation for 
the two-particle Green's function and the associated collective modes. The numerical solution of the 
BS equation in the limit J — >■ shows the existence of two rotonlike collective modes with different 
low-energy Goldstone dispersions and different positions of the rotonlike minima. The two rotonlike 
modes lie outside of the region determined by the lower boundary of the particle-hole continuum, 
and therefore, the two modes are not damped and they should be experimentally observable. In the 
presence of superfluid flow at a certain critical flow momentum, the minimum of the first rotonlike 
mode reaches zero energy, but this occurs before the minimum of the second mode and the lower 
boundary of the particle-hole continuum do, i.e. there are two critical flow momenta related to the 
existence of two rotonlike excitations. 
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I. INTRODUCTION 
A. The Hamiltonian 

One of the key questions of condensed matter physics 
is to understand the nature of single particle and col- 
lective excitations in high-temperature superconductors. 
It is widely accepted that the Hubbard model with an 
attractive on-site interaction plays an important role 
in the qualitative understanding of s-wave supercon- 
ductors, but there is no consensus, however, about 
the theoretical models needed to understand the high- 
temperature superconductivity. The t — U—J model 
with a repulsive on-site interaction can serve as a possi- 
ble model to study high-temperature superconductivity. 
This model was first used in connection with gossamer 
superconductivityjii^ but it also describes the opening of 
a d-wave pairing gap in cuprite compounds, and is consis- 
tent with the fact that the basic pairing mechanism arises 
from the antiferromagnetic (AF) exchange correlations. 
On the experimental side one encounters the problem of 
changing the system parameters t, U and J for a given 
superconductor, because different parameter sets require 
different types of superconducting material. 

Ultracold atomic Fermi gases loaded in optical lat- 
tices present a new opportunity to overcome this obsta- 
cle, and to emulate high-temperature superconductors. 
Optical lattices realize the Hubbard model, if the lat- 
tice potential is sufficiently deep such that the tight- 
binding approximation is valid. In such a system the 
atom-atom interaction can be manipulated in a control- 
lable way by changing the scattering length from the 
Bardeen-Cooper-Schrieffer (BCS) side (negative values) 
to the Bose-Einstein condensation (BEG) side (positive 
values), reaching very large values close to the Fesh- 



bach resonance. On the BEG side of the resonance the 
pseudospin-up and pseudospin-down atoms can form di- 
atomic molecules, and these bosonic molecules can un- 
dergo BEG at a low enough temperature.'^ 

In what follows, we shall focus our attention on the 
BGS side where the existence of a superfluid phase of 
Fermi atoms is expected analogous to superconductiv- 
ity. We shall examine the spectrum of the collective ex- 
citations of population-balanced atomic Fermi gases of 
two hyperfine states with contact interaction, loaded in 
optical lattices. The two hyperfine states are described 
by pseudospins. There are M atoms distributed along 
N sites, and the corresponding filling factor / — M/N 
is assumed to be smaller than unity. For a sufficiently 
deep lattice potential, the system is well described by 
the single-band attractive Hubbard model. 

To make our theory applicable for both superconduc- 
tivity and superfiuidity, the attractive Hubbard model 
will be treated as a J ^ limit of the more general 
t ~ U — J model defined by the following Hamiltonian: 



(1) 



<h3> 



where tij is the single-electron hopping integral, fj, is the 
chemical potential, and rii^a- = V"! o-^i.o- is the density op- 
erator on site i. The Fermi operator ipj ^ {'4'i,cr) creates 
(destroys) a fermion on the lattice site i with pseudospin 
projection a =t,i- The symbol J2<ij> means sum 
over nearest-neighbor sites of the two-dimensional lat- 
tice. The spin operator is defined by ~^ i = (S'f , , Sf) = 
''^i o-'^'^'^'^^o''/2i where "o^ is a vector formed by the 
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Pauli spin matrices {(Jx,ay,az)- The AF interaction 
can be written as JJ2<ij> ~^ j = JJ2<ij> ^i^] + 
5-^E<,,,> [StS- + 5-5+] , where 5± = 5f ± ^Sf . 

B. Collective excitations in moving optical lattices 
in the GRPA 

For the case when the periodic array of microtraps is 
generated by counter propagating laser beams with dif- 
fering frequencies the optical lattice potential is moving 
with a velocity — v (in the laboratory frame) with mag- 
nitude proportional to the relative frequency detuning of 
the two laser beams. In a frame fixed with respect to the 
lattice potential, the fermion atoms flow with a constant 
quasimomentum p — mv, where m is the mass of the 
loaded atoms. For population-balanced Fermi gases the 
order parameter field, ^j{u) — —\U\ < iljj^i{u)ipj^f{u) > 

(or ^*{u) = -\U\ < >): in the mean- 

field approximation varies as $j cx Aexp [2ip.rj]. Here, 
the symbol <> means ensemble average, and A is a real 
quantity which depends on the lattice velocity.— In a 
moving lattice, the formation of BCS superfluidity is pos- 
sible; but due to the presence of quasimomentum p, the 
superflow can break down. 

The stability of balanced superfluid Fermi gases loaded 
into a moving optical lattice has been recently stud- 
ied using the second-order time-dependent perturbation 
theory,- and the Green's function formalism.^ It was 
pointed out that the superfluid state could be destabi- 
lized at a critical flow momentum via two different mech- 
anisms: depairing (pair-breaking) at p^^ and Landau in- 
stabilities at p^^. The depairing takes place when the 
single fermionic excitations are broken, while the Landau 
instability is related to the rotonlike structure of the spec- 
trum of the collective excitations. The superfluid state 
becomes unstable when the energy of the rotonlike min- 
imum reaches zero at a given quasimomentum. The nu- 
merical solution of the number, gap and collective-mode 
equations shows that at a zero temperature the Landau 
instability appears before the depairing mechanismj^ 



C. The Hubbard-Stratonovich transformation and 
the Bethe-Salpeter equation 

It is known that the single-particle excitations of the 
Hamiltonian ([T]) manifest themselves as poles of the 
single-particle Green's function, G; while the two-particle 
(collective) excitations could be related to the poles of 
the two-particle Green's function, K. The poles of these 
Green's functions are defined by the solutions of the 
Schwinger-Dyson (SD) equation G^^ = G'*'^^^ — E, and 
the Bethe-Salpeter (BS) equation [K'^"^-^ - /]* = 0. 
Here, G*^*^^ is the free single-particle propagator, E is 
the electron self-energy, / is the BS kernel, and the two- 
particle free propagator K^^'' = GG is a product of two 



fully dressed single-particle Green's functions. Since the 
electron self-energy depends on the two-particle Green's 
function, the positions of the poles have to be obtained 
by solving the SD and BS equations self-consistently. 

It is widely accepted that the general random phase ap- 
proximation (GRPA) is a good approximation for the col- 
lective excitations in a weak-coupling regime, and there- 
fore, it can be used to separate the solutions of the SD 
and the BS equations. In this approximation, the single- 
particle excitations are replaced with those obtained by 
diagonalizing the Hartree-Fock (HF) Hamiltonian; while 
the collective modes are obtained by solving the BS equa- 
tion in which the single-particle Green's functions are 
calculated in HF approximation, and the BS kernel is 
obtained by summing ladder and bubble diagrams. 

Generally speaking, there exist two different GRPA 
that can be used to calculate the spectrum of the collec- 
tive excitations of the Hubbard Hamiltonian in a station- 
ary (or moving) optical lattice. The first approach uses 
the Green's function method)^"— while the second one is 
based on the Anderson-Rickayzen equations 4ii2r2^ 

The Green's function approach has been used to ob- 
tain the collective excitations in the problem of the ex- 
citon BEC^^* and in s-wave layered superconductors)^ 
According to the Green's function method, the collec- 
tive modes manifest themselves as poles of both the two- 
particle Green's function, K, and the density and spin 
response functions. The two response functions can be 
expressed in terms of K, but it is very common to obtain 
the poles of the density response function by following the 
Baym and Kadanoff formalism, which uses functional 
derivatives of the density with respect to the external 
fields 

The second method that can be used to obtain the col- 
lective excitation spectrum of the Hubbard Hamiltonian 
starts from the Anderson-Rickayzen equations, which in 
the GRPA can be reduced to a set of three coupled 
equations and the collective-mode spectrum is obtained 
by solving a 3 x 3 secular determinant. Recently;^ the 
Belkhir and Randeria 3x3 secular determinant^^ has 
been generalized to the case of a moving optical lattice. 

Instead of mean-field decoupling the quartic interac- 
tion terms, we apply the idea that we can transform 
these quartic terms into quadratic form by making the 
Hubbard-Stratonovich transformation for the electron 
operators. The attempt to decouple the quartic U and 
J terms by using two-component Nambu field operators 

4>{x) = ^jl^l j and ip{y) = (?/'t(y)V'i(y)) does not 

work because the existence of 5+ parts in the J interac- 
tion requires the use of more complicated four-component 

fermion field operators ( ip and ^ obey anticommutation 
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(2) 



In contrast to the previous approaches, such that af- 
ter performing the Hubbard-Stratonovich transformation 
the fermion degrees of freedom are integrated out; we de- 
couple the quartic problem by introducing a model sys- 
tem which consists of a four-component boson field Aa (z) 
(a = 1, 2, 3, 4) interacting with fermion fields ©. 

There are three advantages of keeping both the fermion 
and the boson degrees of freedom. First, the approxi- 
mation that is used to decouple the self-consistent rela- 
tion between the electron self-energy and the two-particle 
Green's function automatically leads to conserving ap- 
proximations because it relies on the fact that the BS 
kernel can be written as functional derivatives of the 
Fock and the Hartree self-energy I = Id + lexc = 
Si:^/SG + S^"/6G = S^<S>/SGSG. As shown by Baymii, 
any self-energy approximation is conserving whenever: 
(i) the self-energy can be written as the derivative of a 
functional $[6], i.e. S = S<i>[G]/SG, and (ii) the SD 
equation for G needs to be solved fully self-consistently 
for this form of the self-energy. Second, the collective 
excitations of the Hubbard model can be calculated in 
two different ways: as poles of the fermion Green's func- 
tion, K, and as poles of the boson Green's function, 
D] or equivalently, as poles of the density and spin 
parts of the general response function, H. Here, the 
boson Green's function, D, is defined by the equation 
D = £)(o)-M)(o)H£)(°) where is the free boson propa- 
gator. Third, the action which describes the interactions 
in the Hubbard model is similar to the action tp'^Aip in 
quantum electrodynamics. This allows us to apply the 
powerful field-theoretical methods, such as the method of 
Legendre transforms, to derive the SD and BS equations, 
as well as the vertex equation for the vertex function, F, 
and the Dyson equation for the boson Green's function, 
D. 

The basic assumption in our BS formalism is that the 
bound states of two fermions in the optical lattice at zero 
temperature are described by the BS wave functions (BS 
amplitudes). The BS amplitude determines the proba- 
bility amplitude to find the first electron at the site i at 
the moment ti and the second electron at the site j at 
the moment t2- The BS amplitude depends on the rela- 
tive internal time ti — t2 and on the center-of-mass time 
(ti+i2)/2>l^ 

exp{z[Q.(r, + r,)/2-tj(Q)(ti +t2)/2]} (3) 
*^.„i(i-»-r,,ti-t2), 




- ^r., p+Q 

k-p,!>-J!'+ ('I'S. „.(p:Si')) 




k -S! P 

k+Q,S!+co'{Q) n'+^Q) 




FIG. 1; Diagrammatic representation of the Bethe-Salpeter 
equation (|4]) for the BS amplitude ^'J?, (k; f2). The single- 
particle Green's function Gn2ni (k, Q.) is denoted by two solid 
lines, oriented in the direction of electron propagation. The 
dashed lines represent the free boson propagator D'^^p{ijj, Q). 

At each of the bare vertices T^a\ni,n2) (represented by cir- 
cles) the energy and momentum are conserved. 



where Q and w(Q) are the collective-mode momentum 
and the corresponding dispersion, respectively. Since 
{rii,n2} = {1,2,3,4}, we have to take into account the 
existence of sixteen BS amplitudes. The Fourier trans- 
form of vI/Q„^(r,0 = / / ^eHk-"-f^*)vI;Q „^ (k; 17) 
satisfies the following BS equation, presented diagram- 
matically in FIG. 1: 



*?.,n,(k;f^) 

= G„,„3 (k + Q,l] + c.(Q))G„,„,(k,17) j — j ^ 

"3 ^5 'k-p,f7-f7'V^. 
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(4) 



Here, G„^„2 (k, fi) is the single-particle Green's function, 
and Id and I^xc are the direct and exchange parts of the 
BS kernel. 

In a similar problem of interacting photons and elec- 
trons in electrodynamics, the direct part of the BS kernel 
does depend on frequency; therefore the solution of the 
BS equation is more complicated. As we shall see in the 
next Section, in the case of the t~U — J model the boson 
propagator ^^^^(k) (a,/? = 1,2,3,4) is frequency inde- 
pendent, and therefore, the following BS equation for the 
equal-time BS amplitude „^(k) = / f vpQ ^^^^ (k; r!) 
takes place: 

^?..u (k) = j ^G„,„3 (k + Q, r! + c.(Q)) G„,„, (k, n) 



{2^Y [-ri°H"3,n5)i?S(k - p)F('')(n6, n,) 
+ ^r(")(n3,n4)i?i'2(Q)F('')(n6,n5)]vI/Q,„^(p). 



(5) 
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In the case oft— U — J model the interaction in the direct 
kernel can be factorized, i.e. of^lk. — p) = /„(k)(7^(p); 
therefore it is possible to obtain the collective excitation 
spectrum using an 80 x 80 secular determinant. Such an 
ambitious task will be left as a subject of future research. 
Instead, we shall discuss the limit J — > in which the 
dimensions of the secular determinant becomes 16 x 16. 

The paper is organized as follows. In the next sec- 
tion we apply the functional-integral formalism to derive 
equations for the single-particle excitations and for the 
two-particle collective modes. In Section III, we numer- 
ically solve the BS equation in the case of ID and 2D 
stationary and moving lattices. We use the mean-field 
approximation for the single-particle Green's functions, 
while the spectrum of the collective excitations is ob- 
tained in the GRPA. The technical details have been de- 
ferred to the Appendix. 

II. FIELD-THEORETICAL APPROACH TO 
THE COLLECTIVE MODES OF THE t~U- J 
MODEL 

A. The functional-integral approach 

The Green's functions in the field-theoretical approach 
are defined by means of the so-called generating func- 
tional with sources for the boson and fermion fields. In 
our problem the corresponding functional integrals can- 
not be evaluated exactly because the interaction part 
of the Hamiltonian ([T]) is quartic in the Grassmann 
fermion fields. However, we can transform the quar- 
tic terms to a quadratic form by introducing a model 
system which consists of a four-component boson field 

Aa{z) {a = 1, 2, 3, 4) interacting with fermion fields ^(y) 



and ipi^)- The spin-dependent nature of the interactions 
requires four spin degrees of freedom of the Bose field 
Aa{z). The action of our model system is assumed to be 

of the following form S = S'^'^ + 5^^^ + S^''-^\ where 
(throughout this paper we use the summation-integration 
convention: that repeated variables are summed up or in- 
tegrated over): 

S'^,'''^ ^\A^{z)Df-\z,z')Ap{z'), 

S^^-^^ ^%)f^^\y,x\ z)^{x)Ao.{z). 

The action S^f^ describes the fermion (electron) part of 
the system. The inverse Green's function of free electrons 
G^^^~^{y\ x) is 4 X 4 diagonal matrix: 

GW-Hy;x) = 

^ exp [l'k.{v, - V,>) - UJrniu - u')] G'(°)-/(k, lUJm), 

kj^TTi f I ( ) 1 

where G\i (k, iWm) = ^^33 (— k, — iWm) — i<^m — 
^(k), and -6*^2^"^ (-k,-ia;,„) = G^"-*"^ (k, itj„) = iujm + 
^(k). The symbol J2uj used to denote /3~^ (for 
fermion fields = (27r//3)(m -)- 1/2); m = 0, ±1, ±2, ...), 
and ^(k) = 2tx{l — cosfcj.) + 2ty{l — cosky) — ^ is the 
non- interacting dispersion on a square lattice,. 

(A) 

The action Sq describes the boson field. The bare 

(A) 

boson propagator in Sq provides the spin-dependent 
interactions U and J, and it is defined as: 



D^°\z,z') = Siv-v') 



USjJ' 



/O 1 
10 


Vo 



j' ,j+a 



( \J- ij \ 

ij |j 

J 

\ J / 



Here the summation on a runs over the nearest- neighbor sites of site j. The Fourier transform of the boson propagator 
is given by 



5(0) (z,z') = l^^e{'[k-(r.-v)— p(---')]}5(0)(k). 



5(°)(k) 



N 

k ujp 

J(k) U-J{\C) 

C/-J(k) J(k) 

2J(k) 

2J(k) 



(6) 



where the symbol is used to denote /3~^ (for boson fields LOp = {2ti / (3)p\p ~ 0, ±1, ±2, ...), and in the case of 
a two-dimensional square lattice J(k) = J (cos kx + cos ky ). 

The interaction between the fermion and the boson fields is described by the action S*'"^'"^). The bare vertex 
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rL (yi;x2 I z) = rL (ii,ui;i2,M2 I = ^(^1 - U2)S(ui - v)5i^i^5i^jT^'^\a) is a 4 x 4 matrix, where 



(7) 



The Dirac matrix 70 and the 4x4 matrices Si are defined as ( <7i are the Pauli matrices): 

^ 



70 



/ 1 
10 

0-10 

Vo -1 



CTyCTiay 



,i = x,y,z. 



The relation between the t — U — J model and our model system can be demonstrated by applying the Hubbard- 
Stratonovich transformation for the electron operators: 



(8) 



/ ^i[A] exp (y)fi") (y; 
= exp [-^%)fi"\y;x\z)${x)D^°^^{z,z')W)^? 
The functional measure -D/i[v4] is chosen to be: 

^i[A] = i?Ae-^-4°W<^"(^^^')^'^(^'),y'^[^] = 1. 

The Hubbard- Stratonovich transformation allows us to map the t ^ U ~ J model onto the model system described 
by the action S. This transformation creates an extra term that can be included in the chemical potential Jl = 
fi— fU/2, and therefore, the mean-field expression for the chemical potential is recovered by the Hubbard-Stratonovich 
transformation. Since there exists a one-to-one correspondence between the t ~ U — J model and the model system 
defined by the action S*, we can obtain the spectrum of the single-particle excitations, as well as the spectrum of the 
collective modes by analyzing the single-particle and two-particle excitations of the model system. 

According to the field-theoretical approach, the expectation value of a general operator 0{u) can be expressed as 
a functional integral over the boson field A and the Grassmann fermion fields ?/; and tl): 



< Tu{0{u)) >= 



1 



Z[J,M] 



DpL[ip,ijj,A]0{u)cyiY> Ja{z)Aa{z) - M{i),'\jj) 



(9) 



where the symbol < ... > means that the thermodynamic average is made, and is an m— ordering operator. The 
functional Z[J, A/] is defined by 



Z[J,M]= D^i[%l),iP,A]0{u)eyi^ Ja{z)Aa{z) - M{i;,ij) 



(10) 



where the functional measure ip, A] = DADipDip exp (S) satisfies the condition J D^lip, ip,A] = l. The quantity 

Ja{z) is the source of the boson field, while the sources Mij{y;x) of the fermion fields are included in the M{ip,'ip) 
term : 

M§, ^) = ^{y)Mn{y; x)^P^ix) + ^pl{y)M2i{y; x)i^t{x) + V'j(2/)Mi2(y; x)^i{x) + V'J (2/)M22(y; x)i;^{x) 
+il;^{y)M3i{y;x)iljf{x) + ipi{y)M4i{y; x)ipf{x) + i/jf{y)M32{y; x)il;i{x) + ipi{y)M42{y; x)il;i{x) 
+ij\{y)Mi3{y;x)x)i;l{x) + ^Pl(y)M23iy;x)^Pl{x) + ,pl{y)^^^^ 

+My)M33{y; x)ijl{x) + ij^{y)M43{y; x)ijl{x) + xjj^{y)M3A{y; x)^|(x) + iJi{y)MM{y; x)^J(x) (11) 

Here, we have introduced complex indices 1 = {ni,xi}, and 2 = {712,2/2} where, xi — {rij,wi}, j/2 = {ri2;^^2} and 
{ni, ^2} = {1,2, 3, 4}. We shall now use a functional derivative 5/6M{2; 1); depending on the spin degrees of freedom, 
there are sixteen possible derivatives. 

By means of the definition (|9]), one can express all Green's functions related to system under consideration in terms 
of the functional derivatives with respect to the corresponding sources of the generating functional of the connected 
Green's functions W[J, M] — In Z[J,M]. Thus, we define the following Green's and vertex functions which will be 
used to analyze the collective modes of our model: 
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Boson Green's function Dap{z,z'): 
This function is a 4 x 4 matrix defined as Dafsiz, z') = jj~^JUJ(P) ■ 

Single-electron Green's function G„j„2 (xi; 2/2): 
This function is a 4 x 4 matrix whose elements are Gn^nj (xi; 2/2) = ^5W/5Mn^n-^^{y2]Xi): 

G(xi;j/2) = 

( <fu U^{xi)^\{y2) \ > <fu f^t(a^i)^I(y2)l > < r„ (7At(a^i)^t(y2)) > < r« (v^t(a^i)^i(2/2)) > ^ 

< fu U^{xi)^{v2)\ > < Tn Ui{xi)xl:l{y2)\ > < (tA^Ix 1)^^(2/2)) > < (V^^(xi)V'i(2/2)) > 

< fu U\{x^)4'\{y2)\ > < fu h\{xi)^l{y2)\ > <fa fvj(a:i)^t(y2)') > < f « f Vj (^^1 (2/2 > 
\^ < f„ (^|(xi)Vj(y2)) > < r„ (^|(xi)V|(y2)) > < Tn (Vl(a:i)^t(y2)) > < f„ (v4(xi)^^(2/2)) > j 



(12) 



Depending on the two spin degrees of freedom, ci and a2, there exist eight "normal" Green functions and eight 
"anomalous" Green's functions. We introduce Fourier transforms of the "normal" Go-j .(jj (k, ui — U2) = — < 

r„(V'tJi,k("i)V'a2,k("2)) >, and "anomalous" Fo-^^cralk, W1-M2) = - < T„(?/'o-i,k(wi)V'o-2,-k(it2)) > one-particle Green's 
functions, where {(Ti,(T2} =tj4-- Here V'^j^('^)7 '0t,k('^) and V'4^j^(u), V'4„k(u) are the creation-annihilation Heisenberg 
operators. The Fourier transform of the single-particle Green's function is given by 



k w„ 



~G(— k, — 



(13) 



Here, G and f' arc 2x2 matrices whose elements are Ga^^a^ and F(j^^„^, respectively. 
The two-particle Green's function K ( '*^' 

\ "2, J/2 ri4,X4 

For different 71.1,71,2,71.3,71,4, there are 256 two-particle Green's functions: 



K 



ni,xi 713,2/3 

712,7/2 714,2:4 



= if 



1 3 

2 4 



(2:1; 2/2) 

(5Af„2„i(?/2;a;i)5M„3„^(y3;a;4) (JM„3„J?/3; 0:4) ' 



(14) 



This definition of K allows us to conclude that if the approximation used for G is chosen in accordance with the 
recipes proposed by Baym and Kadanoff^ii then K is automatically conserving. 

The vertex function rQ(2; 1 | z): 
For fixed a the vertex function is a 4 x 4 matrix whose elements are: 



(^2 



SG„^„^{i2,u2;ii,ui) , , 



SMz') 



(15) 



B. Equations of the boson and fermion Green's functions 



It is well-known that the electron self-energy (electron mass operator) I](l;2) can be defined by means of the 
so-called SD equations. They can be derived using the fact that the measure I?/i[V', t/;, A] is invariant under the 
translations ip ip + Sij] and A-^ A + 5 A: 



D^;^l\z,z')Rp{z') + \Tr (g(1; 2)f(o)(2; 1 | z)) + J„(z) 



0, 



(16) 



G^^l; 2) - G(")-1(1; 2) + E(l; 2) + M(l; 2) = 0, 



(17) 



where Ra{z) ~ 5W / 5Ja{z) is the average boson field. The electron self-energy E is a 4 x 4 matrix which can be 
written as a sum of Hartree and Fock parts. The Hartree part is a diagonal matrix whose elements are: 



S]^(ii,ui;i2,'!i2)„i„2 = if(^°'(zi,wi;i2,W2|j>)„i„2£'^°j(j>;/,uOr^°^(j3,'"3;«4,W4|j',w') 



Gri 



n3"4^ri4n3 



, (74, W4; 73,7^3) 



(18) 
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The Fock part of the electron self-energy is given by: 

„,^ r.(o)/„- „,. „/^^ 



> n4,l4,U4 716,16, We / " > I ! y 



(19) 



The Fock part of the electron self-energy depends on the two-particle Green's function K, and therefore, the SD 
equations and the BS equation for K have to be solved self-consistently. 

Our approach to the t — U ~ J model allows us to obtain exact equations of the Green's functions by using the 
field-theoretical technique. We now wish to return to our statement that the Green's functions are the thermodynamic 
average of the Tu-ordered products of field operators. The standard procedure for calculating the Green's functions 
is to apply Wick's theorem. This enables us to evaluate the T„-ordered products of field operators as a perturbation 
expansion involving only wholly contracted field operators. These expansions can be summed formally to yield 
different equations of Green's functions. The main disadvantage of this procedure is that the validity of the equations 
must be verified diagram by diagram. For this reason we will use the method of Legendre transforms of the generating 
functional for connected Green functions. By applying the same steps as in Ref. Q we obtain the BS equation of 
the two-particle Green's function, the Dyson equation of the boson Green's function, and the vertex equation: 

j^-l f "2,«2,W2 "-3,«3,U3 A ^ ^(0)-l "2,«2,W2 "-3 , «3 , W3 \ _ ^ / ?^2,^2,'"2 n3,i3,M3\ 

n4,i4,U4 J \^ni,ii,ui 77,4,14,7/4 y \^77i,7i,7ti 774,74,774 y' 

D^0{z,z') = D^^l{z,z')+D(^^{z,z'')U^s{z'\z''')Df^{z,z'), (21) 



(22) 



Here, 



p /• ■ I \ r(o)f ■ I ^ . t( n2,i2,U2 773,73,773 \ 

ra(«2,w2;«i,ui I ^)«2ni = r^H«2,u2;ji,ui I 2)„,„i +/ 

Y ''1, '1, "1 "-4, '4, "4 J 

^(0) ( ns, i3, us 776, is, "6 \ p ^ I ^) 

1^774,74,774 775,75,775 J *" ^' ^' '^"""^ 

,^(0) f "-2,«2,U2 n3,«3,'"3\ n r ■ ( ■ ■ \ 

-ft:^ M . . = G„2n3(«2,U3;«2,W2)G„4„i(74,7/4;«l,Wl) 

y"-!,'!,"-! 'I4, '4, "4 / 

is the two-particle free propagator constructed from a pair of fully dressed single-particle Green's functions. The 
kernel / — 6Y./5G of the BS equation can be expressed as a functional derivative of the electron self-energy S. Since 
S — + S^, the BS kernel / = lexc + /d is a sum of functional derivatives of the Hartree Hh and Fock Yip 
contributions to the self-energy: 

712, «2, "2 ns,i'i,Us\ _ ST." {i2,U2]h,U3)n2ni j^f n2,i2,U2 n3,73,773 \ ^ (5S^(72, 772; 73: "3)n2Tti ^23) 



771, 71, 771 774,74,7/4 y (5G„3„^(i3,W3; j4,W4) \^ 771, il, 7/1 774,^4,^4 y (5G„3„, (^3 , "3 ; «4 , "4) 

The general response function 11 in the Dyson equation (PT|) is defined as 

TT / i\ 1^(0) /■ ■ I \ ( n2,i2,U2 773,73,773X^(0).. . I 

Ilaf![Z,z} =r)^'{ll,Ui;i2,U2\ Z)mn2K [ ■ F V (73 , 7/3, l4 , W4 | ^ )»3n4 ■ (24) 

Y '1, "l 'I4, '4, "4 / ^ 

The functions D, K and F are related by the identity: 



1^(0) f n2,i2,U2 ri3,«3,U3 \ p ^- ■ I /N n / ' \ 

^ [ n,,^,,u^ n,,^,,u, ) rM^4, "4; 73, 7/3 | z U^,D0^{z,z) 

( 772,72,U2 n.3,«3,'«3 \ f^(0).. . , is , , x 

= ^ 771, 71, ui 774, *4, u,yV ^^4' "4' "3 u u^.^y^ {z,z), 



(25) 



By introducing the boson proper self-energy P^^iz, z') — n^^(z, z') + Z?|j*j](z, z') one can rewrite the Dyson equation 
ll for the boson Green's function as: 

D-l{z,z')^Dy,~\z,z')-P^p{z,z'). (26) 
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The proper self-energy and the vertex function T are related by the following equation: 



Pais{z,z') = ^Tr f'^°'>{yi,X2\z)G{x2,y3)'^i3{y3,X4\z')G{x4,yi) 



(27) 



It also possible to express the proper self-energy in terms of the two-particle Green's function K which satisfies the BS 
equation K^^ = K'-^^^^ — Ij^, but its kernel Id — 6Y.^ /5G includes only diagrams that represent the direct interactions: 



(28) 



The proper boson self-energy P, the vertex function F and the two-particle Green's function K have common poles. 
Let uJiQ and Q denote the energy and momentum of one of these common poles. Close to ujiq one can write: 



K 



n2,i2,U2 n4,ii,U4 



(29) 



ICJp — U!lQ 

where the amplitudes ip^^ „^ (r^^ , r^^ ; U2 — ui) have the following form: 

V'!5,ni(i"i2:rii;'^2 - Ml) = exp [«Q. (r,, -f J /2] i/^i*^,,,^ (r,^ -r,,;u2 -ui). 
Due to the definition of the bare vertex F*-"-', we have to take into account only the equal "time" u\ — U2 amplitudes: 

^\?,^n, in. - r., ; 0) = ^ ^ exp{zk. in, - r, J}^,|,Q „^ (k), 



where (k) is the equal "time" two-particle wave functions in k-representation. By means of (|29|) and (|28]) we 

obtain: 



where 



EliQ) - JI, and 



lO *IQ 



(30) 



ii,'^^ (0; 0) - 0)) + b^2 (0; 0) - ^4^(0; o) 



bo,z (V'i'llO; 0) - V'4'3(0; 0)) + (5a4 (^/'2'i(0; 0) - "034(0; O)) • The Green's function ^ is invariant under the exchange 
4 o 1 and 2 o 3, therefore, ip^ = -^p^ , and p^ - - -^ t'^— — +1,.+ R- - R-^ r^ - p.. t-i.;o 



means that the proper self-energy has a maximum of six different elements: 

Pii(Q;w) Fi2(Q;w) Pi3(Q;^) Pi4(Q;w) 

p rn-,,^- ( ^i2(Q;^) ^ii(Q;w) -Pi4(Q;c.) Pi3(Q;c<^) 

I Pi3(Q;c^) Pi4(Q;c^) P33(Q;^) P34(Q;c^) 

Pi4(Q;c^) Pi3(Q;c^) P34(Q;w) P33(Q;c^) 



(/?4 . Thus, we obtain that Pn — P22, and P33 ~ P44. This 





b c 






a d 




HI 


d f 


i) 




c h 





(31) 



C. Spin-singlet order parameter 

As we have already mentioned, the BS equation and the SD equations have to be solved self-consistently. The 
so-called P)(°)f(°^ approximation allows us to decouple the above-mentioned equations and to obtain a linearized 
integral equation for the Fock term. To apply this approximation we first use Eq. ([^5]) to rewrite the Fock term as 

S (^1, '^^ij ^2; "^2)711712 — Tq n^l' ^1' ^3; ^3|ji ^)nin3 Dap{j, V; j', w')G„3„4(«3, U3; ^4, 'U4)F^(i4, U4; i2,U2\j', v')n^ri2 , (32) 
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and after that we replace D and V in (|32)) by the free boson propagator D^^^ and by the bare vertex r(°\ respectively. 
In this approximation the Fock term assumes the form: 

E^(il,Ui;i2,M2)nin2 = -ri°H«l ; "1 ; «3 , ^3 1 ; j^ ■"0^^°^ («4, "4; ^2 , "2 |/ 

(33) 

As in the case of the Hartree self-energy, the contributions to mi; 12, M2) due to the elements on the major 

diagonal of the above matrices will be included into the chemical potential Jl. To obtain analytical results, we will 
assume two more approximations. First, we neglect 612 = G21 = G34 = 6*43 = 0, and second, we neglect the 
frequency dependence of the Fourier transform of the Fock part of the electron self-energy , i.e. S-'^(k, iLUm) ~ A(k), 
where the order parameter A is a 4 x 4 matrix. It is known that in order to preserve the antisymmetry in the case of a 
spin-singlet Cooper-pairing, the order parameter A(k) must be an even function of k, so the momentum dependence 
of S'^(k) wiU be 




A(k) = 



In this approximation the Fourier transform of the single-particle Green's function assumes the following form: 



(34) 



G(k, lUr, 



1 



(j^™)2-i?2(k) 



m + e(k) 

lLJ,n + ^(k) 



A(k) 



A(k) 
-A(k) 

-A(k) lLU,n - ^(k) 

iw„ - e(k) 



(35) 



where £;(k) = yF(k)TA(k)2 and ^(k) = 2t^{l - cosk^ 
Green's function in Eq. p3p we reobtain the gap equation: 



2ty(l — cosfcy) — /i. Substituting this single-particle 



^(k) = ^ E [-f^ + 3'^(k - q)] ^ tanh 



(36) 



In the case of the s-wave superfluidity we have J — 0, 
A(k) = A and U < 0, and by using the above single- 
particle Green's function we obtain the gap equation for 
A and the particle number equation for the filling factor: 



/ 



N 
\U 



J2[u\k)fiE{k))+v'ik)fi-E{k))] 



^ f i-E{k)) - f {Ejk)) 
N ^ 2Ep{k) 

1-1 



(37) 



f{x) — [exp (/3a;) -f 1] ^ is the Fermi distribution 
function, and the coherence factors are u^(k) = 
[l + e(k)/£;(k)]/2,«2(k) = l-u2(k). 

The gap equation and the particle number equation in 
the case of a moving lattice are given in the Appendix. 



D. Collective-mode equations 

It is easy to see that in the case of a spin-singlet pair- 
ing a — b = h and c = d = / = 0. From the Dyson 
equation ((26)) , it follows that the spectrum of the collec- 
tive excitations w(Q) could be obtained by setting the 
4x4 determinant det\Df^^^ {Q; uj) — Pq^(Q; uj)\ = 0, or 



r 



det\Sai3 — Daj{Q;i^)Pyl3iQ;i^)\ =0. Thus, we obtain the 
following exact equations for the collective modes: 



= 1 + [L/ - 2 J(Q)] (a -6), 
= 1- U{a + b), 
= l-2J(Q)(a + 6), 



(38) 
(39) 
(40) 



Since the equations (f2T|) and ([24|) hold, the Fourier trans- 
forms of the general response function (1^^ . the two- 
particle Green's function K, as well as the Fourier trans- 
form of the boson Da/s Green's function should share the 
poles defined hy §SI - &■ 

To separate the spin and charge contributions to 
the poles of the general response function we express 
IIq^ (Q; uj) in terms of the Fourier transform of the proper 
self-energy Pai3{Q;uj): 



^S(Q;^)- 



(41) 



n;i(Q;^) = F-/(Q;u 
By means of (|^T|) we obtain 

n(Q; Lu) = n"(Q; uj) + (Q; w 



where n*^(Q;LLi) and n^^(Q;aj) represent the contribu- 
tions to the general response function due to the spin and 
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charge fluctuations: 



~IlP{Q;io) ni'^(Q;w) 

u'^P{Q;uj) 

ni^^(Q;a;) / 

(42) 



U,{Q;lj) U,{Q-oj) 

n,(Q;cj) n,(Q;Lj) 







(43) 



where 



nW(Q;^) = 
(Q;^) = 

nc(Q;w) = 



(a - b)/2 



l + [U-2J{Cl)] (a-b)' 

a — b 
l-2J(Q)(a-6)' 
(a + 6)/2 



(44) 



1 + C/(a + &) ■ 



The coUective modes originating from spin fluctuations 
manifest themselves as poles of the following spin re- 
sponse functions (spin susceptibilities): 



x!(Q;c^) = 2nW(Q;c.) = 
X^(Q;c.)=ni2)(Q;^) 



! + [[/- 2 J(Q)]K''(Q;a;) 



1 - 2J(Q)k^(Q;w) 



(45) 



where «;''(Q;w) = a — b. The second collective mode, 
defined by Eq. pQ]) . manifests itself as a pole of the 
charge density (charge susceptibility): 



X^(Q;c<^) =2ne(Q;c<.) = 



where k'^(Q; uj) ^ a + b. 



1 + 2V^(Q)k<=(Q;w)^ 



(46) 



E. Random phase approximation 

The collective modes can be obtained from Eqs. ((38)) 
- PO)) by solving the BS equation for the two-particle 
Green's function K . The solutions of the BS equation 
are not needed in the RPA, in which K ~ K^^^^ ; or, equiv- 
alently, in Eq. ([27)) for the proper self-energy we replace 
r by f(o): 



RPA 

al3 



-Tr 



f(0)(y,x|z)G(x,y')rf (y',x'|z')G(x',2/) 



(47) 




(b) 



qdlll 



FIG. 2: Excitation spectra in ID stationary optical lattice. 
We set tx ~ t, f — 0.5, U = 2t. The superfluid gap and the 
chemical potential are A = 0.409t and /i = 0.624f. The solid 
line in (a) represents the spectrum obtained by solving the 
equation F\{u},q) — 0, while the dashed line is plotted using 
F2{uj, q) = 0. The spectrum (b) represents the solution of the 
equation _f3(w,g) = 0. 



The RPA expression for the proper self-energy provides 
a-b = h = Ki°''(Q; iLOp) = K'xliQ] lUJp) = Ky°-*(Q; lujp) = 
Kg(Q; lOJp), a + b — Kq{Q; itOp) and c = d ~ f — 0, where 

the general susceptibility k|°-'(Q; zwp), i, j — x,y,z was 



introduced decades ago 

k1?(Q;«Wp) = 



20 



^ ^ ^ Tr |aiG(k; zw„)a_,G(k -I- Q; lujp -f iw„)| 



(48) 



The noninteracting charge and spin susceptibilities 
Kg(Q;a;) and Kg(Q;a;) are: 



Kg(q,iWp) 
Kg(q,iWp) 



-E 



7LU(k,q)[/(^^(k))-/(i?(k + q))] 



k,q 



(*c.,)2-C2(k,q) 



, IV / 1 £(k,q)[l~/(i;(k))-/(i;;(k-fq))] 

N^\mi^\ (zc.p)2-£2(k,q) 

(49) 

Here, C(k, q) = q) - £;(k), e(k, q) = E{V + q) -f 

i?(k), and the following form factors have been used: 

7k, q = "kWk+q + VkWk+q, ^k,q = "k^k+q - Wk^k+q, 
7k,q = Uk?^k-|-q - Mk-l-qWk, ™k,q = "k^k-l-q + "k-l-q^'k- 

(50) 



III. COLLECTIVE MODES OF STATIONARY 
AND MOVING OPTICAL LATTICES AT ZERO 
TEMPERATURE 

In the case of deep optical lattices we neglect the AF 
interaction, and we assume an attractive Hubbard inter- 
action. The collective spectrum is obtained by solving 
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the BS equation at zero temperature in the GRPA. In 
this approximation the corresponding equation for the 
BS amphtude vJ/Q „^ = / g / (Cf*,?.,,.,!^; ^) can be 
obtained from Eq. ([5]) by performing integration over the 
momentum vectors: 



qd/n "'^ 



FIG. 3: Excitation spectra in ID moving optical lattice 
(p=0.21). We set t^ = t, f = 0.5, U = 2t. The superfluid gap 
and the chemical potential are A = 0.420t and /i = 0.655t. 
The solid line represents the spectrum obtained by solving the 
equation F\{u},q) = 0, while the dashed line is plotted using 
F2(a;,(?) = 0. 



71,4 ng 



ni 713 
n2 7^4 



w(Q) X 



f "3 









^ na,n5 J 



(51) 




where the two-particle propagator K^^^^ and the direct 
and exchange interactions are defined as follows: 



"2 ^4 
"2 



/3^/3 '("-4,n-3j 



The boson propagator D^^^ in ([52|) is given by 
limit J(k) = 0. 



(52) 



in the 



FIG. 4: Excitation spectra in 2D optical lattices for stationary 
(a) and moving (b) (px ~ Py = 0.467/\/2) lattices. We set 
= ty = t, f — 0.5, U = 4.5t. The superfluid gap and 
chemical potential in the case of a stationary lattice are A = 
1.334t and /i — 2.234t. For a moving lattice we have A — 
1.386t and ^ — 2.287t. The solid lines in (a) and (b) represent 
the spectrum obtained by solving the equation Fi{uj,q) — 0, 
while the dotted lines are plotted using F2{uj,q) = 0. 



J 



The non-trivial solution of the BS equations ([5T|) exists 
if the secular determinate det\I + UZ\ is zero, where / is 
the unit matrix, and the 16 x 16 matrix Z is given in the 
Appendix. After computing the secular determinant, we 
find that at a given point Q there exist three different 
dispersions ujj{Q,),j = 1,2,3, which are the solutions of 
the three equations Fj{ujj, Q) = 0, where: 



Fl{uj, Q) = 1 + {In + + Iy,-y) U + + Il^ll„ 



J jZ 



7^2 



(53) 
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F2(w,Q) = ! + (-/„ 



2^,7) - 



7^2 



T T'' 



j2 < .T J T 

r2 o T_ T _ r I t 



7,7 ^ 7,7/ 



7,7''7,i 



(54) 



(55) 



The definitions of b and Ja,b are given in the Appendix, 
and we have used the relation /s; = Ii i + Im,m — ^7,7- 

The dispersion of the first rotonlike mode a;i(Q), ob- 
tained by solving numerically the equation _Fi(wi, Q) = 
has been previously obtained by the equation-of-motion 
method)^ and by locating the poles of the density re- 
sponse function!^ The numerical solution of the equation 
F2{(j-'2, Q) = provides the existence of a second roton- 
like collective mode cj2(Q)- The two rotonlike modes 
have different low-energy Goldstone dispersions and dif- 
ferent rotonlike minima. 

In FIG. 2 - FIG. 4 we have plotted the excitation spec- 
trum of superfluid Fermi gases in ID and 2D optical lat- 
tices assuming the same strength of the interactions and 
filling factors as in FIG. 4 (?7 = 2t, / = 0.5) and as in 
FIG. 9 ([/ = 4.5^,/ = 0.5) in Ref. [|. In the presence 
of superfluid flow, as the flow momentum increases, the 
two rotonlike spectra lean toward the left side and the en- 
ergies of the two rotonlike minima decrease. It is worth 
mentioning that the two rotonlike modes lie outside of the 
region determined by the lower boundary of the particle- 
hole continuum; therefore the two modes are not damped 
and constitute propagating modes and they should be ex- 
perimentally observable. At a certain flow momentum, 
the minimum of the first rotonlike mode reaches zero en- 
ergy, but this occurs before both the minimum of the sec- 
ond mode, and the lower boundary of the particle-hole 
continuum do. At this critical flow momentum the first 
mode is destabilized due to the spontaneous emission of 
rotonlike excitations; but to destabilize the second mode 



one has to increase the flow momentum. Thus, there ex- 
ist two critical lattice velocities which are determined by 
the existence of two rotonlike modes. 

The third mode, 0^3 (Q), is obtained by solving the 
equation i^3(a;3,Q) = 0. This mode lies entirely in- 
side the particle-hole continuum Therefore, it is a non- 
propagating mode and cannot be an experimentally ob- 
servable mode. 



IV. CONCLUSION 

In conclusion, we have studied the superfluid state of 
Fermi gases in deep optical lattices. Using the Bethe- 
Salpeter equation in the GRPA we have obtained the col- 
lective mode spectrum of the attractive Hubbard model 
in the presence of superfluid flow. We found that the 
spectrum of the collective excitations has two rotonlike 
modes with different low-energy Goldstone dispersions 
and different rotonlike minima. Both modes are experi- 
mentally observable because they are separated from the 
region determined by the lower boundary of the particle- 
hole continuum. As the flow momentum increases, the 
two rotonlike spectra lean toward the left side. At a 
certain critical flow momentum, the energy minimum of 
the first collective mode hits zero, but this occurs before 
both the minimum of the second mode, and the lower 
boundary of the particle-hole continuum do. 
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Appendix 

The BS equation ((5T|) written in the matrix form is (j + U "i! ^ 0, where / is the 16 x 16 unit matrix. The 
16 X 16 matrix Z is defined as follows: 

/ 



Z = 
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v 

The transposed matrix of ^ is given by: 

i]/^ — \ \I>*^ vt*^ O*^ vJ/*^ vl/^ vl/^ vl/'^ vl/*5 vl/*-^ vl/Q flfr'^ \If*^ \If*^ 

^ ~ 1 *1,2 *1,3 *1,4 *2,1 *2,2 *2,3 *2,4 *3,1 *3,2 *3,3 *3,4 *4,1 *4,2 *4,3 *4,4 



At zero temperature the elements of matrix Z are: 

2 {Im.m + -^7.7 '^J'j^m) j -^4444 — 2 



^1111 = 7^ {Im,m + — '2Jj,m) , -^4444 = 2 (^"i,™ + -^7,7 + '2'Jj,m) , ^1144 = 2 (^'.' + -^7,7 + '^Jl,l) 

K4411 — — {IlJ + Ij^j — '2'J-fd) , Klin = 2 [^h,m + -^7,7 + ^/,7 ^ '^'.m) , -^^^4414 = 2 (^'.'" ^ -^7.7 '^':7 ~ "A,™) 
Kiiii — - {IlJ ~ Ij^j) , A'i444 — — {Il^m + ^7,7 + Jl,^ + J'f,m) , A1411 — - { — Il^m — ^7,7 + Jl,^ — J^,m) , 



where the following symbols are used 

Ia,b 



J_ Qk,Q&k,Qg(k, Q) 
w2_e2(k,Q) ' 



JaJ 



-E 



w^ak,Q&k,Q 



Here, a and & are one of the form factors (|50|) . 

The above relations can be easily extended to the case of a moving lattice. In the presence of quasimomentum p 
the following definitions and relationships hold: 



up(k) = Wi[l + x(k;p)/£;p(k)], «p(k) = Ji[l-x(k;p)/£;p(k)], i?p(k) = Jx'(k; p) + , 



X(k; p) = - K(p + k) + e(k - p)] , 77(k; p) = - K(p + k) - C(k - p)] , i?±(k; p) - 77(k; p) ± i?p(k). 
At zero temperature the gap equation for Ap and the particle number equation are: 

l^ x(k;p) 1 
A £:p(k) ' A ^ 2£;p(k) ' 

In the presence of quasimomentum p the elements of the 16 x 16 determinant Z are the same as in the case of a 
stationary lattice, but the definitions of laj, and J^.h ^ire as follows: 

^ 1 <Q^k.Qgp(k,Q) J 1 <Q^k.Q[^ + ^p(k,Q)] 

A^- [,, + f]p(k,Q)]2-£2(k,Q)' A^[c. + J7p(k,Q)]2-e2(k,Q)' 
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where £p(k, Q) = £'p(k + Q) + EpCk) and np(k, Q) = ?7p(k) — ?7p(k + Q). a and b are one of the following form 
factors: 

7k,Q = Wp(k)Mp(k + Q) + wp(k)t;p(k + Q), /P Q = up{k)up{k + Q) - Wp(k)wp(k + Q), 



7pk, Q Up(k)wp(k + Q) - Wp(k)up(k + Q), q Up(k)up(k + Q) + Up(k)Mp(k + Q) 



r 



It is well-known that the singularity of the integrals la.t 
and Ja,b at energies uj — ep(k, Q) — r2p(k, Q) corresponds 
physically to the possibility of depairing into two fermion 
excitations with energies E+{k+ Q; p) and — £^_(k;p). 
The spectrum for this kind of excitation is known as 



the particle-hole continuum.— The lower boundary of 
the particle-hole continuum is defined by the condition 
mirik (-E+(k + Q; p) — i?-(k; p)) where the minimum is 
to be taken over all the possible values of k. 
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